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Abstract. We construct a new finite difference method for the flow of 
ideal viscous isentropic gas in one spatial dimension. For the continuity 
equation, the method is a standard upwind discretization. For the mo- 
mentum equation, the method is an uncommon upwind discretization, 
where the moment and the velocity are solved on dual grids. Our main 
result is convergence of the method as discretization parameters go to 
zero. Convergence is proved by adapting the mathematical existence 
theory of Lions and Feireisl to the numerical setting. 
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1. Introduction 

In this paper, we will develop a convergent finite difference method for 
the flow of an ideal viscous isentropic gas in one spatial dimension. We will 
assume that the flow may be modeled by the Navier-Stokes system (cf. [16]): 

ft + H = 0, inM+x(0,L), (1.1) 

(gu) t + (gu 2 ) x = fiu xx - p(g) x , in R + x (0, L). (1.2) 

The unknowns in this system are the fluid density g = g(t, x) and the fluid 
velocity u = u(t,x). For an isentropic flow, the ideal pressure law takes the 



Date: March 13, 2013. 

2010 Mathematics Subject Classification. Primary: 35Q30,74S05; Secondary: 65M12. 

Key words and phrases, ideal gas, isentropic, Compressible Navier-Stokes, Conver- 
gence, Compactness, Finite elements, Discontinuous Galerkin, Finite Volume, upwind, 
Weak convergence, Compensated compactness. 

The research is funded by the Research Council of Norway (proj. 205738). This research 
was conducted at Paul Sabatier University, Toulouse. The author wishes to express his 
gratitude towards Paul Sabatier University for its hospitality. 

1 



2 



KARPER 



form 

p = ag 1 , a > 0, 

where the value of 7 is determined by the specific gas in question. In this 
paper, we will be forced require that 

^<7<2, 

to prove convergence of the method. Note that this significantly limits the 
physical applicability of the convergence result. While monoatomic gases 
(7 ~ |), such as helium, are included, diatomic gases (7 ~ |), such as air, 
are not. The condition is necessary for convergence, but not stability, of the 
method. 

At the boundary, (1.1) - (1-2) is augmented with no-slip conditions, 

u(t,0) = u(t,L) = 0. 

and the initial conditions, 

< q(0,x) = q (x) £ L 7 (0,L), u(0,x) = u (x) £ L°°(0,L). (1.3) 

Existence of global classical solutions for the system (1.1)-(1.2) is well- 
known and was first established by KaneP [9] for smooth data (see also [14]). 
For non-smooth initial data, corresponding results have been obtained by 
David Hoff [1, 7]. These works are all based on the same approach: find 
pointwise upper and lower bounds on the density, then use these bounds 
to estimate the second derivative of the velocity. With smooth initial den- 
sity, the density remain H 1 , while with discontinuous initial density, the 
density is at most BV as initial discontinuities persists for all time (see 

]). The pointwise bounds on the density are obtained by tracking certain 
quantities along streamlines rendering the entire existence theory essentially 
Lagrangian. 

In the literature, one can find a huge variety of numerical methods appro- 
priate for (1.1)-(1.2). However, very few of these methods have been proven 
to converge. This lack of rigorous results are most likely a consequence of 
the Lagrangian nature of the existence theory. In fact, prior to this paper, 
all convergent methods have been discretized in Lagrangian coordinates and 
are due to David Hoff and collaborators [18, 19, 20]. That being said, there 
are also several existence results that utilize discrete Lagrangian approxima- 
tion schemes (e.g. [1, 7, 8]) to construct solutions. For practical applications, 
results in Eulerian coordinates is often desirable and a mapping from dis- 
crete Lagrange to discrete Euler is costly. For this reason, most practitioners 
would employ an Eulerian discretization. 

In this paper, we will discretize the equation in Eulerian coordinates. As 
a consequence, obtaining pointwise bounds on the density becomes highly 
involved and will not be pursued in this paper. Instead, we will develop a 
convergence theory in the spirit of the continuous existence theory [4, 15, 17] 
as sparked by P. L. Lions. In particular, we will not obtain any form of 
continuity of the discrete density and instead prove strong convergence of 
the density using renormalization and what is known as the effective viscous 
flux (cf. [17]). 



CONVERGENT FINITE DIFFERENCES FOR ID VISCOUS ISENTROPIC FLOW 3 



In more than one spatial dimension, the literature is almost void of con- 
vergent numerical methods. For the full system, the only result is the paper 
[13], by the author, in which a convergent finite element discretization of 
(1.1)- (1.2) is developed. Over the last years, there have also been developed 
a series of convergent methods [2, 3, 5, 10, 11, 12] for the Stokes version of 
(1.1)-(1.2). The method we shall develop and analyze in this paper can be 
considered as an archtype for all the methods mentioned above. That is, in 
one spatial dimension, all of these methods are in a sense equivalent and of 
the form we shall consider here. In addition, the method we present here is 
the first finite difference method for which convergence is proved. 

The paper is organized as follows: In the next section, we will define the 
numerical method and state the main convergence result. Then, we will 
derive stability of the method and some other basic properties. In Section 
4, we will derive an equation for the effective viscous flux and also provide 
a higher integrability estimate for the density. Section 5 concerns passing 
to the limit in the method. In particular, we prove that the limit almost 
satisfies (1.1)-(1.2). The remaining ingredient is to prove strong convergence 
of the density which is proved in the final section. 

2. The numerical method and main result 

Our method will be posed on a uniform grid in both space and time. In 
time, we shall approximate at discrete points t k = kAt, where At is assumed 
to be of the order Ax. In space, we will divide the domain [0, L] into N 
intervals of length Ax = L/N. It will be convenient to write this as 

[0, L] = |J[xi_i/2, x i+ i /2 ], i = 0, . . . ,N - 1, 

i 

where we have introduced the slightly confusing notation 

Xi-\/2 = iAx, i = 0,...,N. 
We shall also need the dual grid given by the midpoint nodes 

Xi = (i + i j Ax, i = 0,. . . ,N- 1. 

The numerical method will approximate the density on the dual grid and 
the velocity on the standard grid: 

g(kAt,Xi) « q\, i = 0,...,N-l, 

u(kAt, Xj_ 1/2 ) « u-L 1/2 , i = 0, . . . , N. 

Staggered grid of this kind is necessary for incompressible flows, and as a 
consequence also widely used to develop all speed compressible flow methods. 

2.1. Discrete operators. To discretize the convective terms, we shall uti- 
lize an upwind method. To facilitate this, we introduce the notation 

u + = max{u,0}, u~ = min{u, 0}. 

We shall also need the average velocity over an interval 

U,-_1 lo + It; 4-1 /9 
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For the continuity equation (1.1), we shall use the following upwind flux 

Up(£U)i+l/2 = Qiutrl/2 + Qi+lS+1/2- 

This upwind flux will also be used for the momentum equation, where we 
upwind the (averaged) momentum 

\Jp(guu) i+1/2 = (QiUi)uf +1/2 + (g i+1 u i+1 )u7 +1/2 . 

The remaining derivatives will be discretized using the operators 
a p fi+1 -fi a v i+l/2 ~ v i-l/2 

8*11*/ = —^-, 9 iV = , 

which defines the following (standard) Laplace operator 

a a a u i-l/2 ~ 2n i+ i/ 2 + u i+3/2 

^i+l/2 u h — O i+ i/ 2 OiU h — ^2 • 

For time discretization, we will use implicit time-stepping 



£K fK 

ak f J i J i 



2.2. Numerical method. The numerical method is defined as follows: 



Definition 2.1. Given At, Ax < 1, and initial data (1.3), define the 
numerical initial data 

Qi = -*r- / 6o(y) dy, u i _ 1/2 = u (x i _ 1/2 ), i = 0, . . . , N - 1. 

X J x i-\/2 

Determine sequentially the numbers 

(f?f,«f), i = 0,...,N-l, k = 0,...,M, 
solving the nonlinear system 

d? 6i + d i \jp(e k u k ) = o (2.1) 
k ( gm + Qi+m+A ( Uy(Q k u k u k ) i+3/2 - \Jp(g k u k u k ) i _ 1/2 \ 
9t V 2 ) + ^ 2A^ ) 

= /A +1/2 4 - d l+1/2 p(g k h ), (2.2) 

U-l/2 = u N-l/2 = °" 



It is not completely trivial that our numerical method is well-defined. 
Since the system of equations (2.1)-(2.2) are both nonlinear and implicit, 
existence of a solutions needs to be established. 

Lemma 2.2. Let < At, Ax < 1 be fixed. The nonlinear implicit system 
of equations (2.1)-(2.2) admits at least one solution. 

Proof. The existence can be proved using a toplogical degree argument. 
Since the proof is very similar to the corresponding result in [6, 11, 13], 
we do not give the details here and instead refer the reader to any of the 
mentioned papers. □ 



From [6] and [1 1] , we also have that the method preserves strict positivity 
of the density: 
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Lemma 2.3. Let {(g k , u k _^ / 2 )}j satisfy the continuity scheme (2.1). Then, 

k ^ • k— 1 f 1 \ 

mm > mm ft ' 1 



1 + At maxj |m* +1 / 2 | 



2.3. Extension. To analyze the finite difference method, it will be of great 
convenience to extend the numerical solution to all of [0, T) x (0, L). Since the 
density appears nonlinearly in the pressure, we will use piecewise constants 
to extend it: 

g h (t, x) = g k , V(t, x) E [t k , t k+1 ) x fo_ 1/2 , x l+1/2 ). (2.3) 

For the velocity, we shall utilize piecewise continuous linears in space and 
piecewise constants in time. 

u h (t,x) = u k + [ X l l ~ 1/2 ) U k - u k ) , (2.4) 



i-l/2 I ^ J V i +1 / 2 ~ i " 1 / 2 

for all (t,x) E [i fc ,i fc+1 ) x [^-1/2) ^1+1/2)- Note that an extension of the av- 
erage velocity u k is now given by the 1? projection onto piecewise constant: 



u h (t,x) = IL%[u h ] = — / u h dx, x E (x i ^ 1 /2,x i+1 / 2 ). 



2.4. Main result. The following theorem is our main result in this paper. 

Theorem 2.4. Assume we are given initial data {qo,uq) satisfying (1.3) 
and a finite final time T > 0. Let {(Qhi u h)}h>o be a family of numerical 
solutions, constructed through Definition 2.1 and (2.3)-(2.4), with 

At = Ax = h. 

Then, as h — > 0, Uh — 1 u in L 2 (0, T; Wq' 2 (0, L)), qu — > g a.e on (0,T) x 
(0,L), where (g,u) is a weak solution of (1. 1) -(1.2) : 

T rL r-L 

I g(4>t + u- 4>x) dxdt = / £o<XO, •) dx, 
io Jo 



/ / (£>u)uf + gu 2 v x + iiu x v x - p(g)v x dxdt = 
Jo Jo Jo 

for all (<t>,v) EC7 °°([0,r) x (0,L)). 



L 



Theorem 2.4 will follow as a consequence of the various results stated and 
proved in the upcoming sections. The proof will be completed in Section 6. 



3. Stability and energy estimates 

In this section, we we will derive a numerical analog of the continuous 
energy estimate. This estimate yields in particular stability of the method, 
but it will also provides us with necessary L p bounds uniform in the dis- 
cretization parameters. 
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3.1. Rernormalized continuity scheme. We shall need the following 
renormalized continuity scheme at several occasions. The term renormal- 
ized is motivated by the corresponding continuous equation and its role in 
the existence theory (cf. [15]). 

Lemma 3.1. Let B G C a (0,i2 + ) and define b(z) = zB'(z) - B(z). If 
{(gf, u^ +1 , 2 }i,k satisfies the continuity scheme (2.1), then the following iden- 
tity holds 

d k t B{ ei ) + di Up (s(eV) + K&diu" + v [s(gj),« fc ] = o, (3.i) 

where V is given by 



V 



B(g\),u = AtB"(g*)\d* 6i 



B"(g\)[g i+ i - gifu. 1/2 + B"(g\)[gi - ^_i] 2 u H 



Here, g*, g\, and g\, are some numbers in the range [g\~ X , g 1 ^ 1 ], [^j^f+i]? 
and [^f_!,^], respectively. 

Proof. We begin by multiplying (2.1) with B'(g) to obtain 

= UpG?u). 

By applying Taylor expansion, we find that 

d?B(g) + AtB" (g*)\dj£ Qi\ 2 = -B'ig^ Vp(gu), (3.2) 

where g* is some number in [^f -1 , gf]. 

By adding and subtracting and applying another Taylor expansion, 

(Ax)B , (g l )d l \Jp(gu) 

= B'(g i )g i (u i+1 / 2 - Wi-1/2) + - ft]w^.i/2 

+ - £>i-i]ut 1/2 

= B'(gi)gi(u i+ i/2 ~ Ui-1/2) + [#(£>i+i) - B (^)] n r+i/2 (3.3) 
+ [-B(ft) - 5(^_i)]it+ /2 + S"(4)[gi+i - ft] 2 « 



+ 

i— 1/2 

= (Ax) ^B'(gi)gidiU + d i+1 / 2 B(g)uJ +1/2 + d i _ l/2 B(g)ur_ l/2 

+ B"(eJ)[ft+i - £>i]\>i/ 2 ~ - £ , i-i] 2 ^_i/ 2 

Here, g\ and g| are some numbers in [£j,£>j+i] and [f?i, f?i-i], respectively. 
Next, we calculate 

B'(gi)gidiU + 5 i+1/2 S(^)u7 +1/2 + ^_i/ 2 S(g)n+ 1/2 

= b(gi)diU + B(gi)diU + d i+1/2 B(Q)u7 +1/2 + di_x/ 2 B(g)ur_ 1/2 (3.4) 

= b(gi)diu + di Up(B(g)u) 

We conclude the proof by combining (3.4) in (3.3) and (3.2). □ 
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3.2. The convection operator. To prove stability of the method, we shall 
multiply the momentum scheme (2.2) by u k +1 , 2 an d sum over au i- F° r this 
purpose, we shall need the following identity for the convective discretiza- 
tion. 



Lemma 3.2. The following identity holds 



N-l 



i=0 



Up(g fc W) i+3/2 - Up(g fc S fc n fc )^i/ 2 \ k 

2Ax J Ui+l ' 2 

N-l /i-Sfc|2\ 
i=0 ' v 7 



Qi 



7fc|2 



+ No- 



where the numerical diffusion term Mi is given by 

„ AT 1 

(Ax 



M 2 



i=0 



Up ( Q k U k 



i+1/2 



0; 



+1/2U 



Proof. By applying summation by parts, we see that 



iV-l 



AxJ2 



\Jp(e k u k u k ) i+3/2 - Up(g k u k u%_ 1/2 \ k 



i=0 



2Ax 



U i+l/2 



N-l 



-Ax ^v{Q k u k u k ) i+1/2 



i=0 
N-l 



U i+3/2 U i-l/2 
2A~X 



k u k u k ) i+1/2 d i+1/2 u k h . 



i=0 



(3.5) 



Next, we apply the definition of Up(-) and add and subtract to deduce 
\Jp(g k u k u k ) i+1/2 d t+1/2 u k h 



k -~~k 



u i+l) U i' + l/ 2 



di+l/2U h 



Qi u 7+l/2 + ^i+l n i+l/2J + ) di+l/2 u % 



+ ^H+i/2 + ^+i<;i /2 ] («? - ) ^+1/2^ 



i+l/2 



fc„,& 



Up £Tli 



i+1/2 



/2 



+ 1/2 
fe|2 



Ax 
~2~ 



Up ( e V 



i+l/2 



<9i+l/2^ 



We then conclude the proof by setting this identity in (3.5). 



□ 



3.3. Energy estimate. We are now ready to prove stability of the method. 
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Proposition 3.3. Let {(g^ , u^_^ / 2 )}i,fc; be the numerical solution obtained 
through Definition 2.1. The following stability estimate holds, 



max 



^ It;" 1 ! 2 1 

AxY J Qr- J ^ L + — PWi 



7-1 



M 



i=l 



k=l i i- 

= Ax]T^ + -l-^), 



(3.6) 



where the numerical diffusion terms are 

M 

k=l i 
M 



M = (At) 2 Ax^^p"(^)|5 4 V 

k=l i 
M 

Af 2 = At(Ax) 2 *£p"U) d l+1/2 Q h 



k=l 
M 



It 



i+1/2 



Jfc-l 



Af 3 = (AtfAxY J Y,\- ^ 

k=l i 
, 2 M N-l 



Ma 



k=l i=0 



Up ( Q k U k 



i+1/2 



9 i+ l/2U 



Proof. We begin by multiplying (2.2) by u k +1 ^ 2 Ax and sum over all i to 
obtain 



Ax ^<9, 



£W + Qi+Wi+i \ k 

U i+l/2 



-AxJ2 



Vp(g k u k u k ) i+3/2 - Vp(Q k u k u k ) l _ 1/2 \ k 



2Ax 



U; 



i+1/2 



AxJ2\ 9 i+l/2 u h ~P(0i)di 



where we have also applied summation by parts together with the boundary 



condition u 



1/2 ~ a N-l/2 



0. Next, we apply Lemma 3.2 and Lemma 3.1 



(with B{z) = t^t[P(z)) to deduce 



k ( Qi u i + Qi+l u i+l \ k 



U, 



-Af 2 - AxJ2 \di+i/2Uh 



i+1/2 u t 
2 1 



d k t Qi 



7fe|2 



(3.7) 



7-1 



$p{Qi)-vW)M 



To proceed, we observe the following identities 



Ax^d, 



QiUi + Qi + lU i+ l \ k 



U ■ 



i+1/2 



Ax^<9 t fc (QiUi) '' 



(3.8) 
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and 



dt (QiUi) u 1 ! - d k tQi 



Tfc|2 



dt I ft 



U; 



(3.9) 



+ 



2 ) 2At 
By applying (3.8)-(3.9) in (3.7), we deduce 



u k -u k ~ l 



Ax^ 



k io m 
t \ 6i 



+ 



+ la, 



'i+l/2Uh 



2At 

2 



u k -u k ~ l 



+ M2+V[p(g), 



+ -^d k p( 0i 
7-1 

= 0. 



We conclude by multiplying with At, summing over k = 1,...,M, and 
recalling the definition of V[-]. □ 

The previous stability estimate also provides some uniform integrability 
estimates on various quantities. To state these, let us introduce the notation 

fh e b lp(Q), 

to denote the case when fh is bounded in LP uniformly with respect to both 
At and Ax. 

Corollary 3.4. Let (gh,Uh) be the numerical solution constructed through 
Definition 2.1 and (2.3)-(2.4). Then, 

g h G 6 L°°(0,T;i7(0,L)), p(g h ) G& L°°(0, T; L 1 (0, L)), 

u h G L 2 (0, T; <' 2 (0, L)) n L 2 (0, T; L°°(0, L)). 
As a consequence, 

Q h u h € b L°°(0,T;L^(0,L)), 6h \u h \ 2 G& L°°(0, T; L 1 (0, L)), 

g h u h G fe L 2 (0,T;L 7 (0,L), Qh \u h \ 2 G b L 2 (0, T; (0, L)). 

Proof. From Proposition 3.3, we have that Qh Gb ^°°(0, T; L 7 (0, L)), Uh Gb 
L 2 (0,T;Wq' 2 (0,L)), and g h u 2 h G b L°°(0, T; L 1 (0, L)). The remaining bounds 
follows directly from these using Sobolev embedding and interpolation esti- 
mates. □ 

3.4. Control in time. We end this section by establishing some weak con- 
trol on the time-derivatives. For this purpose, let 

f h (-) - f h (- - At) f k - f- 1 



for all Jfc = 1,...,M. 



At 



At 



t G [kAt, (k + l)At), 



Lemma 3.5. Let (gh, Uh) be the numerical solution constructed through Def- 
inition 2.1 and (2.3) -(2.4). Then, 



d h t Qh Gi L 2 (0,T;W-^(0,L)) , 
d[ l (g h u h ) Gi L 1 (0,T;W-^(0,L)). 



(3.10) 
(3.11) 
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Proof. 1. We first prove (3.10). Let <fi S C°°(0, L) be arbitrary and define 

I f x i+l/2 



Ax L 



4>{y)dy, i = 0,...,N-l. 



-1/2 



Now, multiply the continuity scheme (2.1) with Axfa, sum over all i, and 
apply summation by parts, to deduce 



AxJ2(dh 



Ax VUp(A fe ) 1+ i/2 . v 

< c||^IIl^(o,l)II' u /iIIl°°(o,l) 



i+1 — <Pi 



(3.12) 



L^T(0,L) 

The last inequality is an application of the Holder inequality. Now, since 
(3.12) holds for all <f>, we can conclude that it continues to hold for all 

1 7 (0 L) 

(j) £ W 't _1 ' By multiplying with At and summing over all k, we obtain 



< l|£ , h||L° o (0,T;L7(0,L))l|«/i||L 2 (0,T;L°°(0,L))) 



£ 2 (0,lW-i.7(0,.L)) 

which is (3.10). 

2. Let € Wq'°°(0, L) be arbitrary and define 

^-1/2 = ^(^-1/2), i = 0,...,iV. 

By multiplying the momentum scheme (2.2) with Axf i+1 / 2 , summing over 
i, and applying summation by parts, we obtain 



Ax ^2 9 t(.Q$i + Qi+iUi+i)v i+l/2 

i 

= Ax^Up(£ fe S fc u fc ) 

i 

+ A x J2»dr kfVi+1/2 ~ Vi - 1/2 



k~k„.k\ ( V i+3/2 ~ Vj-1/2 

2Ax 



Ax 



k\ f v W/2 - v i-l/2 

Ax 



>(0,L) ll U *ll L ^ (0 , L) 



+ ^\\( u h)x\\L2(0,L) \\Vx\\ L 2 {0jL) + ||p(^)||l-(0,L) KIIl-(0,L) 

Next, we multiply with At and sum over all k to conclude 



At Yl 



Ax^dt{QiUi + Qi + lU i+ l)v i+l / 2 



< C\\v x \\l°°(0,L) (ll n ^llL 2 (0,T;L°°(0,L))ll^lli° o (0,T;L 7 (0,L)) 

\\( u h)x\\L 2 {0,T;L 2 (0,L)) + \\p(.Qh)\\L°°(0,T;L 1 {0,L))) ^ ^11^ IU°°(0,L) i 

where the last inequality is Corollary 3.4. Finally, we calculate 



dt(QhUh)v dx 



dt 
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S^ahf ~n r i+1/2 f v i+i/2 + v i-i/2 , ,\ , 
2^ d t / = - v{x) dx 



'i+l/2 



+ Ai£ 



2^ 9 i / - v ( x ) dx 



< C (l + ||£/i%||l1(0,T;L1(0,.L))) II u 2;||l°°(0,L); 
form which (3.11) follows. 



□ 



4. The effective viscous flux 

The purpose of this section is to derive an equation for the quantity 

F h = n(u h ) x -p(g h )- 

This quantity is often termed the effective viscous flux and stands at the 
center of both the available existence results (in more than ID) [15, 4] and 
the available numerical results (cf. [13]). Specifically, both higher (than L l ) 
integrability on the pressure and strong convergence of the density is proved 
using the upcoming equation (see Proposition 4.2 for details). 

To derive the desired equation, we shall need the following discrete Neu- 
mann Laplace operator 

-Aiq h = -did i+ i /2 qh = fi, 
qn = Qn-i, 
9-1 = gl- 
and we observe that — Aj is nothing but the standard 3-point Laplacian on 
the density grid. Moreover, the Neumann condition is realized by adding 
a shadow cell on both sides of the domain (0,L). Due to the Neumann 
condition, Aj is only well-defined for sources fh of zero mean. 

In the upcoming analysis, we shall not need Aj directly, but its discrete 
derivative: 

fh = -di+i/2lh> 

where we observe that the Neumann condition renders 0jV-i/2^~ [/*] = 
and cL^A" 1 [/^] = 0. We shall also need the inverse of the A i+1 / 2 operator 
occurring in the momentum scheme (2.2): 

- A i+1/ 2 M = W i+l/2, 

where Wh solves the linear Dirichlet system 

-Aj+i/2^/1 = -d i+1/2 diw h = v i+1/2 , w_ 1/2 = w N _ 1/2 = 0. 

We will mostly be interested in the discrete derivative 

diA7^ 1/2 [v h ] = -diw h , 
The following result follows from standard summation by parts. 
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Lemma 4.1. The following duality holds, 

Aa? E^+Va^+i^Aj" 1 [h] = -^xJ2diA7^ 1/2 [v h ]fi. 

i i 

Proof. By direct calculation, using the definition of Aj and Aj +1 / 2 , and the 
respective boundary conditions, we deduce 

Ax^2v i+1 / 2 d i+1 / 2 Aj l [fh] 

i 

= a^^^a^h) (di+t/zK 1 [f h ]) 

i 

= -Ax^2(diA-^ /2 [v h ]) {d l d t+1/2 A- 1 [f h ]) 

i 

= -AxJ2(diA^ 1/2 [v h ])f i , 



which concludes the proof. 



□ 



The following proposition gives the effective viscous flux equation we shall 
need in the convergence analysis. The two error terms appearing in (4.1) 
will be bounded below. 

Proposition 4.2. Let (gh, Uh) be the numerical solution constructed through 
Definition 2.1 and (2.3) -(2.4). For any m = 1, . . . , M , we have that 



k=0 i 



AtAx^^2\J V { e k u k u k ) 



i+1/2 



k=0 i 



M 
L 



-A*J>A7 + i /2 

i 

+ A*J>A7 + 1 1/2 



of 1 



0-0 , -o 
Vi u i + Si+l u i+l 



} + E? + E™, 



(4.1) 



where the numerical error terms E™ and E™ are given by (4.6) and (4.8), 
respectively. 

Proof. Let m = 1, . . . , M be arbitrary. By multiplying the momentum 
scheme (2.2) with v i+l / 2 (AtAx), where 



v i+i/2 = d i+l/2Al 



1 


" k M~ 

& ~ T_ 


, M= [ 






Jo 



g° dx, 
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and summing over all i and k = 0, . . . , m, we obtain the starting point 



k=0 i 



1/2 



k=0 i 



+ At Ax ^2 E 



QjUj + Qj + lU i+ i 

2 j V i+l/2 

'Vp{g k u k u k ) l+3/2 - \Jp(g k u k u k ) i _ 1/2 \ k 

V i+l/2 



k=0 i 



2Ax 



5i + S 2 . 



(4.2) 



Here, we have introduced the quantities 



Sl = AtA X £ E # ( 9m + v k +1/2 

k=0 i ^ ' 



fc=0 i 



2Ax 



n k 

V i+l/2 



We will now rewrite the sums Si and S 2 using the definition of w f +1 / 2 an< ^ 
the continuity scheme. We will treat the least involved term, namely S 2 , 
first. For this purpose, we apply summation by parts and the definition of 
v i+i/2 t° calculate 



S 2 = -AtAxY,Y. V ^S k u k u k ) t+1/2 

k=0 i 
in 

= -AtAx^^Up(/u fc n fc ) i+1/2 



J i+3/2 ~ v i-l/2 

2Ax" 



(4.3) 



k=0 i 
in 



+ AtAsEE U Pte fe « fe A+i/2 . 

k=0 i ^ ' 

For the S\ term, we first apply summation by parts in time followed by 
an application of Lemma 4.1 to deduce 



k=0 i 



h>i u i ' fj+i 



,S, AtArY.Y. ' 

k=0 i 

Til 

k=0 i \ 
m 

= AiA*£5>A-\ /2 
AxE^ A m/ 2 



a i+ i/ 2 A- 1 



L 



fc-i^fc-i , fc-i^fc-i 

Bi u i + Qi+i u i+i 



(4.4) 



Bi a i + Bi+i u i+i 



+ AxE^A" 1 



Bi a i ' Bi+i a i+i 



b¥ 



+1/2 



Bi 
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Next, we multiply the continuity scheme (2.1) by AtAxqf, where 



^7+1/2 



"err 1 



^ Ui+l u i+l 



and sum over all i and k = 0, . . . , m to obtain 



1/2 



fc=0 i 



fc-i s fc-i , k-lqk-i 



+1/2 



A:=0 i 
fit 



fc-l~fc-l fc-l-fc-1 
"i ^ Ui+l u i+l 



AtA X E E U P ( AV1/2 g - + 

fe=0 t V / 

AtAx 53 E u P ( AVi/i. * 4 2 t+1 m + 



(4.5) 



where the last identity follows from summation by parts, the definition of 
A 7, 1 ,_, and where we have introduced the error term 

1+1/2' 



-AtAxY^YUp( k u k ) l+1/2 



k=0 i 



k-l^k-1 , Jc-l-fc-l 



(4.6) 



Qi u t + g i+1 



To conclude the proof, we shall need to further rewrite (4.5). By adding 
and subtracting, we obtain the identity 



up(A fc ) m/2 



QM + Q k +iu k i+l 



(k + k k — k 

Si U i+l/2 + Si+l U i+l/2 



' uM + Ui+iUj+i ^j 



n k v + ' k v k ( ^ + g * +1 I I n k v + ' k n k ( ^ +1 
Qi u i+l/2 u i 2 I BiU i+1/2 Q i+1 I 



+ n fc v~' k v k ( 6 * + ^ +1 
"t- ft + in i+1 ^ 2 'Uj +1 i 2 



&+l w i+l/2& 2 



+ QiQi+i 



7/ fe 

U i+l/2 
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Consequently, by combining this identity (4.5) and (4.4), we conclude 

Sl = AtA*f;£Up(,W) ( e ^fA 

k=0 i V / 

i 

+ A *E 5 A r + i/ 2 



-1 

i+1/2 



6i U{ -r Q i+ iU i+1 



0-0 , cfl 



qI + e^ + e?, 



where 



i+1/2 



E™ = AtAxYY.£e k i+1 

fc=0 i 

Finally, we apply (4.3) and (4.7) to (4.2) to discover 

m 

^eeKi^j - 9 i+i/2P(Qh) 

k=0 i 

= AtA*f;]rup(/sV) m/2 (^) 



V i+l/2 



n mcrra , m cm 
Qi U i T Qi + \U i+l 



+ Ax^9,A, H 1 1/2 

i 

This and a final summation by parts, 



0^0 , syo 



^EEK+i^ - d i+i/2P(Q k h) 



k=0 i 



-At Ax Y E i^ 9 ' 



k=0 i 



J i+l/2 

M 
L 



concludes the proof. 



(4.7) 



(4.8) 



□ 



4.1. Bound on the error terms in (4.1). In order for the previous propo- 
sition to be useful, we will need to prove suitable bounds on the error terms 
E\ and E\ ■ It is at this stage we will need to impose some requirements 
on the value of 7. No other results in the paper imposes any unphysical 
restrictions on 7. 

We begin by deriving a bound on E\ . 

Lemma 4.3. Assume that the adiabatic exponent satisfies 

4 



7 > 



3' 



and let E^ be given by (4.8). There is a constant C > 0, independent of At 
and Ax, such that for any m = 1, . . . , M , 

\E?\ < (Ax)^, 
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Proof. By two applications of the Cauchy-Schwartz inequality and a stan- 
dard inverse estimate, we deduce 



H+l/2 



|£ 2 H<|£f|<AiAxEE^+i 

k i 

< (Ax)C ( AtAx E E \ diU i+i 



u i+i ~ u i 



/2 



k i 



x H u /ilU2(0,T;i«>(0,L))l|e/i|lL«'(0 ) T;i 4 (0 I X)) 



37-4 



< (Ax) 1 (Ax) V3 T ; C< (Ax)^C : 



where the norms are bounded due to Corollary 3.4. 



□ 



We are now ready to bound the other error term in (4.1). The following 
lemma is the sole reason for the requirement 7 > |. 

Lemma 4.4. Assume that the adiabatic coefficient satisfies 



7> 2' 



and let be given by (4.6). There exists a constant C > 0, independent of 
discretization parameters, such that for any m = 0, . . . , M, 



27-3 



\E?\ < h~C. 



Proof. Prom (4.6), we have that 



m 

\2a™\^ Y^ 



(AtyAx^^Vp( 6 k u k ) l+1/2 

k=0 i 



dt(6iUi) + d£(Q i+1 u i+ i] 



(4.9) 



Let us now examine one of the terms in E\. By adding and subtracting, we 
write 



(At) 2 Ax J2 E Up( A*W#(ft%) 

k=0 i 
m 

< (At) 2 A*EE| U P^ n *W^ \ Q tQi 



k=0 i 



(4.10) 



+ (At) 2 AxEE| u P^ V W2 e k i~ l ^ 



k=0 i 



^:S 1 + S 2 . 
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To bound the Si term, we multiply and divide by y/p"(g\) and apply the 
Cauchy-Schwartz inequality 



Si = (At) 2 Ax EE |Up(A*W6?| \%Qi\ 

k=0 i 

/ m 2 ^ 

< (At)*C AxA*EE |Up(A*)i+i/2^ { 8 \f^ 
V fc=0 i / 

/ m 2 \ 2 

x (a*) 2 ax£J>"(^) # ft 



(4.11) 



fc=0 j 



Next, we make several applications of the Holder inequality together with 
standard inverse estimates, to deduce 



Ax At E | U P( A fc )m/2£- 



k=0 i 



<T II II 2 II 1 1 2 — -y II i|2 

^ \\ U h\\L°°(0,T;L 2 (0,L))WSh\\ L°°(o,T;L°° (0,L)) H H (0,T;L 2 (0,L)) 

2-7 2— J 

< C(Ax) -' ||0fc||£oo( O ,T;2,T(O,i)) 

x (Ax) 2 (5-^)||^n h || 2 ^ 

L°°(0,T;L7+I(o,L)) 

-,,27-3 

< (Ax) _1+ ~C7, 



(4.12) 



where we have used Corollary 3.4 to conclude the last inequality. 
Combining (4.12)-(4.11) and recalling that At = Ax, yields 



27-3 

Si < (As) "ST. 



(4.13) 



Next, we turn to the S2 term in (4.10). An application of the Cauchy- 
Schwartz inequality yields 



S 2 = (At) 2 Ax J2 E H e kyk 

k=0 i 
/ m 

< (At) 2 AxEE^-i 



Jfe-1 



k=0 i 



(4.14) 



(m 
AtAxEE| U P(^Wi 



/2 



k=0 i 



< (At) 2 C|| QhUh\\L2(0,T;L2{0,L)) \\Qh 



L°°(0,T;L°°(0,L)) 
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Next, we proceed as in (4.12) to discover 



AtAx^Yl PP(A*)i+i/2 



k=0 i 

^ II l|2 || || 

S \\QhUh\\L2^ 0!T . L 2^ 0tL ^\\Qh\\L°°(0,T;L°°(0,L)) 



2(I_i) 2 _l 

< C(Ax) V2 -<J\\g h u h \\ L 2 {0tT . LJ ^ L)) (Ax) T\\Qh\\L°°(0,T>L-y(0,L)) 

1,27-3 

< (A X y 1+ —c, 

(4.15) 

where we have used Corollary 3.4 in the last inequality. 

By combining (4.14) and(4.15), and recalling that At = Ax, we conclude 

S 2 < (Axf^TC. (4.16) 
By setting (4.13) and (4.16) in (4.10) we obtain 

Si + S 2 < (Ax)^C, 

and hence we have the desired bound for the first term in (4.9). The second 
term in (4.9) can be bounded by the exact same arguments. □ 



4.2. Higher integrability on the density. From Corollary 3.4, we only 
know that p(gh) is uniformly (in At and Ax) bounded in L°°(0, T; L 1 (0, L)). 
Hence, it is unclear whether p(Qh) actually converges to an integrable func- 
tion. In the following lemma, we prove that the pressure has more integra- 
bility than provided by the energy estimate. 

Lemma 4.5. Let (Qh,Uh) be the numerical solution constructed using Defi- 
nition 2.1 and (2.3)-(2.4), with 

3 

7> 2- 

There is a constant C > 0, independent of discretization parameters, such 
that 



[ [ qI +1 dxdt < C. 
Jo Jo 



a 



7 + 1 



In other words, p(g h ) e L°°(0,T;L -r (0,L)). 
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Proof. First, we rewrite the equation (4.1) in the form (4.1) 
fT r L 

p{Qh)Qh dxdt 



J 



7 

o Jo 



M f T f L 

p(Qh)~j- dx + fi / (u h ) x Q h dxdt 

L Jo Jo 

m 

+ At Ax Yl ^v{Q k u k u k ) i+l/2 

k=0 i 

- A ^E^ A r+i/2 



® 



n Mc.M , M cjM 



1/2 



o-o , o 



To bound the terms involving the discrete inverse Laplacian, we shall use 
the elementary bound 



SjAr, 1 



i+l/2 



Sh u h 



Together with the Holder inequality and Lemmas 4.4 and 4.3, this readily 
provides the bound 

rT r -L 

p(Qh)Qh dxdt 



j 



< All, 

^ £ H6/Hl£,oo( 0l T;£T(0. 



n 



£,)) + - IK n ft)a:llL 2 (0,T;L 2 (0,L)) 



17+1 



+ e || Qh IIl7+1(o,T;L^+ 1 (0,L)) + £ P hn ^ 1 1 L 1 (O.TjL 1 (0,L)) 

2 7 -3 

+ A^C||^%|| i o O(0)T . i i (0iL)) + (Ax) 2 ^ C. 
The proof is completed by fixing e sufficiently small. 



□ 



5. Weak convergence 

In this section, we will pass to the limit in the numerical method and 
prove that the limit is almost a weak solution to the compressible Navier- 
Stokes equations. Our starting point is that Corollary 3.4 allow us to assert 
the existence of functions 

u G L 2 (0, T; Wq' 2 (0, L)), 

£eL°°(0,T;LT(0,L)), ( 5 .1) 

^£1^(0,^^(0,1)) 
and a subsequence hj — > 0, such that 

Qh^gm L°°(0,r;LT(0,L)), 

u h ^umL 2 {0,T;W£' 2 {0,L)), (5.2) 
pW^^ini^(0,r;L^(0,i)). 
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Note that we cannot make the identification p(g) = p(g) as this would 
require the density to converge strongly. We will prove that this is indeed 
true in the next section. 

To conclude convergence of the product terms appearing in the method, 
we shall need the following lemma from [11]: 

Lemma 5.1. Given T > and a small number h > 0, write [0, T) = 
U^ijtfc-i.tfc) with t k = hk and Mh = T. Let {fh}™ >0 , {9h}^ >0 be two se- 
quences such that the mappings t \— > gh(t, x) and t t-t fh(t, x) are constant on 
each interval (tk-i,tk] and assume that {fh}h>o, {gh}h>o conver 9 es weakly 
to f andg in L? 1 (0,T; L qi (0, L)) and L p ' 2 (0, T; L 92 (0, L)) , respectively, where 
1< Pl , gi < oo and i + i = i + i = 1. // d\g h G b ^(0, T; W- 1 ' 1 ^, L)) 
and h a \f(-,x + h) - f(-,x)\ £5 L P2 (0,T; L q2 (0, L)), for some a < 1, then 
9hfh ~ 1 gf i> n the sense of distributions on (0, T) x Q. 

Lemma 5.2. Given the convergences (5.2), 

QhUh QU in L 2 (0, T; L 7 (0, L)), 

g h u h ^ gu in L°°{0,T;L^{0,L)), (5.3) 

QhUhUh, Qh\uh\ 2 , Q h \u h \ 2 -± qu 2 m L 1 (0,T;L^(0,L)). 

Proof. From Corollary 3.4 and Lemma 3.5, we have that 

u h G b L 2 (0,T;W Q 1,2 (0,L)) and d£ eh € b L 2 (0, T; PF _1,7 (0, L)). We can then 

apply Lemma 5.1, with = Qh and fh = Uh, to conclude 

QhUh QU in L 2 (0,T;L 7 (0,L)). 

Next, we notice that 

2 



\uh ~ u h \\ 2 L 2 {0:L) - ^2 







fXi+l/2 




Ax 


'xi-1/2 





Uh(y) dy - u h (x) 



dx 



< (Ax) 2 \\(u h ) x \\ L 2( 0iL) , 

where the last inequality is the Poincare inequality. Hence, by writing 

QhUh = QhUh + {QhUh- u h ), (5.4) 

and passing to the limit, we conclude the second convergence of (5.3). 
From Corollary 3.4 and Lemma 3.5, we have that 

Qh\u h \ 2 Gb L 1 (0,T;L^(0,L))and^(^) G b L X {Q, T; W -1,1 (0, L)). Lemma 
5.1 is then applicable, with gh = QhUh and fh = Uh, yielding 

Q h u h u h QU 2 in L l (0, T; Lr+i(0, L)). 
The remaining convergences can be obtained similarly (i.e (5.4)). □ 

5.1. Convergence of the density scheme. We now prove that the limit 
(g, u) is a weak solution of the continuity equation. 

Lemma 5.3. Let (Qh,Uh) be the numerical approximation constructed using 
Definition 2.1 and (2.3)-(2.4). The limit (q,u) is a weak solution of the 
continuity equation. That is, 

Qt + (qu) x = 0, 

in the sense of distributions on [0,T) x [0,L]. 
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Proof. 1. We first claim that (gh,Uh) satisfies the equation 

I [ dt(Q h )<P-Q h u h cf> x dxdt = P 1 (<j ) ), (5.5) 
Jo Jo 

for all <p G Cg°([0,T) x [0,L]), where P x (4>) = At£ fc P-f(0) is given by 
(5.10). 

To prove the claim (5.5), let <p G C^°([0,T) x [0,L]) be arbitrary and 
define 

I r(k+l)At rx i+1/2 

^ = a7^/ / 4>{t,x) dxdt, i = 0,...,N-l. 

AtAx J kAt J Xt _ 1/2 

Now, multiply (2.1) with ip^Ax, and sum over all i, to discover 

Ax Y^(. d t Qi)4>i = Ax £ fe u P (A fc : 



= £ Up(/u fc ) m/2 (0* - 0(A*fe, x i+1/2 )j (5.6) 

- Up(A fc )i-i/ 2 - ^(AtMi-1/2)) • 
To proceed, we add and subtract to derive the identities 

Up(A*0i+i/2 (V- - ^(AiMi+1/2)) 

= eN*+i/2 (fi ~ H&tk, x i+1/2 )) (5.7) 
+ (Qi+l ~ &i) u 7+i/2 (&* ~ 0( A * fe > x i+i/2)) . 

and 

Up(A fc )i-i/2 - 0(Atfc,x f _i /2 )) 

= eNti/ 2 - <f>(Atk, Si_ 1/2 )) (5.8) 

- (e? - ?ti)«t 1/2 - <t>(&tk,Xi-i/2)] 

By applying (5.7)-(5.8) in (5.6), we obtain 

Aa?£(#ft)$ = £ ^u-+i /2 - </>(Atk,x i+ 



~ 6i u ti/2 (<f>i ~ 0(A^,^_ 1/2 )) 



"i-1/2 

We then apply Green's theorem to obtain 
Ax£(^)^ 



£ / (W* " dx + PjM), 

i Jxz-l/2 



(5.9) 
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where we have defined 



p£{</>) = Ax^ [(V1/2/) u+_ 1/2 + (d i+1/20 k 

i 

To proceed, we observe the two identities 



Vl/2 



k i 



JO 



(dtQh)4> dxdt 



(5.10) 



(5.11) 



and 



At EE 



T /-L 
JO 



QhUh<t>x dx, 

(5.12) 



fif (jWi -<£(AtA;,x))J rfa; 

/ Qh(uh)x(<Pi ~ 4>( x )) ~ QhUh4>x dx = - 
o Jo 

where the last equality follows by definition of <\>\ 

By combining (5.9)-(5.12), we obtain (5.5), which was our claim. 

2. Next, we prove that the Pi{4>) converges to zero as h — > 0. More 
precisely, we claim that 

\Pl{4>)\<hkc. (5.13) 
To prove this, we apply the Cauchy-Schwartz inequality to obtain 

Fx (<£) = AxAtJ2J2 [(Vi/2£ fe ) «£i/2 + (?i+i/2Q k ) V+l/2 



x -^(Atfc.Xi-i/a)) 
<c(At(Ax) 2 ^^p"(4) aw-i/afig 



I "4+1/2 1 



(5.14) 



(Ax)" 1 /* (AtAx^J> m/2 |(^ 



fc i 



l 

2\ 2 



0*-0(Atfc, 2^1/2) 

Note that the first term after the inequality is bounded by the energy esti- 
mate (3.6). Now, the Holder and Poincare inequalities provides the bound 



(Ax)- 1 / 2 AtAz^EKn/alGif 



fe\2-7 



k i 



<p(Atk, x i _ 1 / 2 ] 



< (Ax) 2 C[| V0[|£, a o( OjT .x«.( Q)i )) [|«^[| 22 (0)T;L oo( 0iL) ) lbfe|U°°(0,T;LT(0,L)) 

< (Ax)5C||V0|| iO o (Oir;Z/ oo (OiL)) , 

where we have used Corollary 3.4 to conclude the last bound. Together with 
(5.14), this proves our claim (5.13). 
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3. Let us now send h — > in (5.5) and thereby conclude the proof. For 
this purpose, we shall need the following elementary identity 



dt(Qh)<fi dxdt = - 

JO JO 



f Q h (t- At)9 t fe dxdt - [ g° h 4>(At) dx, 
o Jo 



where we have also used that (f>(T, •) = 0. With this identity, (5.5) tell us 
that 



— / / Qh(t — At)dt<f) dxdt — QhUh 4>x dxdt 
Jo Jo 

= [ q14>{M) dx + Pi(<f>), e c °°([o,r) x L\X). 

Jo 

From (5.13), we have that Pi((f>) — > as h — > 0. Hence, there is no problem 
with sending h—}Q, using (5.2) and (5.3), to conclude that 

Qt + (qu) x = 0, 

in the sense of distributions. This concludes the proof. □ 



5.2. Weak limit of the momentum scheme. In this subsection, we pass 
to the limit in the momentum scheme to conclude that (g,u) is almost a 
weak solution of the momentum equation. We begin by deriving an integral 
formulation of the momentum scheme (2.2). 

Lemma 5.4. Let (g^, u^j be the numerical solution constructed through Def- 
inition 2.1 and (2.3)-(2.4). Then, for all sufficiently smooth v, 



t rL 

dt(6hUh)v - gh(u h ) 2 v x dxdt 



o Jo 

h) ./ii 

where the numerical error term is given by 



T r L 

ip(Qh) - K u h)x)v x dxdt + P 2 (v) 



(5.15) 



" !+1/2 ~J V i~l/2 + V i+l/2 



k i Jx i-1/* V J 

At 

+ T 



i 

~ ^+1/2(^+1/2-^-1/2) • 

Proof. Let v G Cq°([0, T) x 0, L) be arbitrary and introduce the notation 



«tl/2 = V 



(kAt, Xi _ 1/2 ) , i = 0,...,N, k = l,...,M. 
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Now, multiply (2.2) with v k +1 ^ 2 Ax, and sum over all i to obtain 
/\ T 

~Y z2 d t(eiUi + ^+1/2^+1/2)^+1/2 

i 

+ \ £ (up - up »U < 5 ' 16 > 

= Ax J] (/xA i+1/2 u fc - d;+i /2 p(^)) vf+i/2- 



We will now write each integral term in (5.16) in the form (5.15). Let us 
begin with the right-hand side. 

1. Using summation by parts, we calculate 



>i+l/2P\Qh)) V i+l/2 



V i+l/2 ~ V i-l/2 



(5.17) 



AxY, [^A i+1/2 u k -di. 

i 

i 

= (p(0i) - M« fc ) / + ' v x {kAt,x) dx 

L 

{p(Qh) ~ Li(ui) x )v x {kAt,x) dx, 





which is of the form we wanted. 

2. For the time derivative term, we perform summation by parts to 
deduce 



J^dtiQiUi + ft+l/2^+1/2)? 



i+1/2 







dt(g h u h )vdx (5-18) 

\- r i+i/ % fe , ^f-1/2+^+1/2 , A , 

L 

dt(QhU h )v dx + H\ 
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3. For the last term, we begin applying summation by parts 



i 

= -\ ^2^p(e k u k u k )(v i+3/2 - Vi_ 1/2 ) 

i 

= ~\ ^Z^{Q k u k u k ){v i+V2 - v i+1/2 ) 

i 

\^Z^v{Q k y k u k )(v l+1/2 - Ui-i/a). 



(5.19) 



2 

Now, by adding and subtracting, we develop the identity 

Up(Q k U k U k )(v i+3/2 -V i+1/2 ) 

= ^+1^+1^+1/2(^+3/2 - v i+i/2) (5.20) 
- (ft+inf +1 - p 4 fe uf) u+ + 1/2 K +3/2 - u m/2 ). 
Similarly, we derive the identity 

Up(/«*n fe )( %1/2 -Ui-1/2) 

= ^Sfu i+1/2 (^ +1/2 - Uj_i/ 2 ) (5.21) 
+ - u~ +1/2 (v i+1/2 - Vi_ 1/2 ). 

By applying (5.20) and (5.21) in (5.19), we obtain 

\ E (up {^ k - k ) l+3/2 - up (^ w L 1/2 ) * l/a 

i Jx i-l/2 

+ \Y, (^+l 2 tl - ft* [«tl/ 2 ( V <+3/2 " «i+l/ 2 ) (5 ' 22) 
i 

~ u ffl/a( u i+l/2 - Wi-1/2) 

=: ~y q\ (u^y ' v x dx + H k . 

4. By applying (5.17), (5.18), and (5.22), to (5.16), multiplying with At, 
and summing over all fc,we discover 



[ [ k (g h u h )v - Qh (u h ) 2 v x dxdt + At V (e\ + E\ 
Jo Jo 1. 



A: 

I: 

r rL 

{p{Qh) ~ K u h)x)v x dxdt, 



10 Jo 
which is (5.15) with 



P 2 (v) = -AtJ2(Hi+H : 
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This concludes the proof. □ 

Next, we prove that the error term in (5.15) converges to zero as h — > 0. 

Lemma 5.5. Let P2(v), be as in the previous lemma. There is a constant 
C > 0, independent of h, such that 

\P 2 (v)\ < /l 5 C||w a; ||£, o(o i T;i«'(0,L)) 5 
for alive L°°(0,T; W 1,oo (0,L)). 
Proof. By definition, we have that 



p 2 (v) = -Atj2Y,[ 

1, .- J X 



OtiQiUi 



V i-l/2 + V i+l/2 



v(x) dx 



k i Jx *-l/2 



At 



k i 



k~k 



+ T S E {-■ ■ I " ~> "'' 



^1/2(^+3/2 " Wi+1/2) (5.23) 
^+1/2(^+1/2 - ^-1/2) 



=: Si + S 2 . 

Let us now bound the Si and S 2 term separately. 



1. Bound on S\: By adding and subtracting, we write 

^1+1/2 



5i = -a*eE 

fc i 



fe i ■ / *»-i/a 

^1+1/2 



fc ! 



U i-l/2 + ^1+1/2 



1/2 

x i+l/2 
i-1/2 



V i-l/2 + U j+l/2 



' fc 1 k 

V i-1/2 + U i+l/2 



u(x) dx 



v{x) ] 



(5.24) 



-:Ki+K 2 . 



To bound the Aj term, we apply the Cauchy-Schwartz and Poincare in- 
equalities to obtain 



1*1 



x i+l/2 



V i-l/2 + V i+l/2 



v(x) I dx 



fc i 



fc i 



^-1/2 + U i+l/2 



(5.25) 



v(x) dx 



1 

2\ 2 



< (At) 2CA^2(A2;)||t; x || L oo {0iT;L 2 (0iL)) < h2C\\v x \\ Lx ,^ T . L oo^ L)) , 

where we have used that the term involving d\ Uh is bounded by the energy 
estimate (3.6). 
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To bound the K2 term, we apply yet another application of the Cauchy- 
Schwartz and Poincare inequalities to discover 



\K 2 



k i 



x i+l/2 



' k 1 k 

V i-l/2 + V i+l/2 



v{x) dx 



-1/2 



k i 



(5.26) 



^2- 7R fc|2 



U f-l/2 + U i+l/2 



u(x) <ix 



1 

2\ 2 



< (At) 2 C , ||u /l || L 2 (0)T;Lt x ) ( 0iL)) A^2(Ax)||w :c || L oc {0iT;L cx 1(0iL)) . 
Setting (5.25)-(5.26) in (5.24) yields 

Si < /i I C||?; :I .|| L oo( 0iT;i 2 oo ( 0>i )). 



(5.27) 



2. Bound on S2: By adding and subtracting, we calculate 
At 



So 



k i 



n i+l/2(^+l/2 — ^—1/2) 



< (Ax) 11% 11^(0^00) AxAt^^ 1^+1/2^1 



(5.28) 



/2ll U i 



+ |^i+l/2 u k| 1^+1/2 1 



=: (Ax) \\l°°( Q ,T;L°°) ( K 3 + K*) ■ 

To bound the K3 term, we apply the Cauchy-Schwartz inequality, followed 
by the Holder inequality, and obtain 



K 3 = At Ax Yl \ d i+i/28h\ l«i+i/ 2 l 



k i 



< 



(At)"* (At) 2 Ax^^p"(^)|3 m/2 ^ 



2 I.. A 



W 



i+X/2\ 



k i 



x (AtAx^^K^H^I 2 ^) 2 -' 

V k i / 

< (At)-lcxi||^||| 3(0iT;ioo(0iL)) 

< (At)-iC7(At)l(5-i)||^||| 2(0iT;ioo(0)L)) < (At)-ic, 

where we have also utilized a standard inverse estimate in time and the 
energy estimate (3.6). 



(5.29) 



28 



KARPER 



Next, we apply the Holder inequality to deduce 
K 4 = AtAx^Y^ \ 9 i+i/2Uh K+i/2\0i+l 

k i 

< II(^/i)2;||l 2 (0,T;L 2 (0,L))II u /i||l 2 (0,T;L° o (0,L))II^/iIIl°°(0,T;L 2 (0,L)) ( 5 ' 3 °) 

<(Ax) 2 1 \\Qh\\L°°(0,T;L-r{O,L))C, 

where we have utilized an inverse estimate and the energy estimate (3.6). 
By setting (5.29) and (5.30) in (5.28), we conclude 

\S 2 \ < llvrcllioo^o^.^^x,)). (5.31) 
Finally, we apply (5.27) and (5.31) in (5.23) to conclude the proof. 

□ 

Lemma 5.6. Let (Qh,Uh) be the numerical approximation constructed using 
Definition 2.1 and (2.3)-(2.4). The limit (g,u,p(g)) satisfies 

(gu) t + (qu 2 ) x = nu xx - p(g) x , (5.32) 

in the sense of distributions on [0,T) x 0,L. 

Proof. We begin by writing (5.15) in the form 



T i-L 

(QhUh)(-At, •) d^v + Qh{uhf v x dxdt 

o Jo 

T 

(p{Sh) - K u h)x)v x dxdt + / (g° h u° h ) v(A, •) dx + P 2 (v). 

o Jo 



(5.33) 



The convergence (5.2) provides p(gh) p(q) and from Lemma 5.5, we have 
that -P 2 ( u ) — J- as /i —7- 0. Moreover, from Lemma 5.2, we know that g\{dh — i 

qu in L°°(0, T; L 2j j + 1(0, L)) and g h \u h \ 2 gu 2 in L x (0, T; (0, L)). 
Hence, there is no problems with sending h — > in (5.33) to conclude the 
proof. □ 

5.3. The effective viscous flux limit. We are going to end this section by 
passing to the limit in the effective viscous flux equation (4.1). To achieve 
this, we shall need the following result. 

Lemma 5.7. Let {fh}h>o, {9h}h>o be two sequences satisfying 

• For each fixed h > 0, fh and gh are piecewise linear and piecewise 
constant, respectively, with respect to our grid with At = Ax = h: 

t . X — Xj^ 1/2 

9h(t,x) = g,- L , v h (t,x) = v { _ l/2 H — {v i+1/2 - %-i/ 2 ), 

for (t, x) e [kAt, (k + I) At) x (^-1/2,^+1/2), VA, i. 

• Ash->0, A A / andg h ^ g in L°°(0, T; L p (0, L)) and 
L°°(0,T; (0, L)), respectively, where l/p+ 1/q < 1. 

• The discrete time derivate of fh satisfies 

d*v h e h L x {Q,T;W- l > x (S),L)). 
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Then, for any t G (0,T), 



lim AiV^A" 1 



h^O 



+1/2 



d x A- l [f(t)]g(t) dx, 



(5.34) 



where k is given by k = [t/h\ and Ad denotes the Laplace operator with 
homogenous Dirichlet conditions. 



Proof. Let be the extension of diA^^,^ 
q h (t,x) = diAT^ 1/2 



V i+l/2 



to all of [0,T) x (0,L); 



V i+l/2 



V(t, x) G [&Ai, (/c + l)Ai) x (a;i_i/2,^i+i/2) an d relevant i,k. From the 
requirements on fh, it is clear that 

Then, we can apply Lemma 5.1, with fh = 9h = Qh to conclude that q\ — i 
in the sense of distributions. Thus, % — > q a.e as h — > 0, where the 
convergence may take place along a subsequence. Observe that the limit 
must satisfy 

q = d x A D 1 [v], 

where A^ is the Laplacian with Dirichlet conditions. 
Next, let qf^ be the linear- in-time interpolation of q^: 



Qh(t,x)=q k h (x) + 



t - kAt ( k+1 



At 



q* h +i (x)-q k h (x)), te [kAt, (k + 1) At), 



for k = 0, . . . , M — 1. Since we have that 



d -^^-=dUh Gb L\Q,T;L\0,L)), 

we have that qfc — > q in C(0, T; L p (0, L)) and sup t \\qfc — qh\\LP(o,L) ~~ 0) f° r 
any p < oo. Finally, we write 



/2 



and send /i — > to discover (5.34). 



Qh9h dx+ i (q% - q h )g h dx, 



□ 



Equipped with the previous lemma, we now pass to the limit in the effec- 
tive viscous flux equation (4.1). 

Lemma 5.8. Assume that 7 > § and that the convergences (5.2) and (5.3) 
holds. Then, for any t £ (0, T), 



Jo 



inu x -p(g)) ( g 

t r L 

qu 2 dxdt 



M 
L 



dxdt 



o 



10 JO 

where the overline denotes the weak L 1 limit. 



L 

d x A~j^ [gu] g dxdt 



(5.35) 



s=0 
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Proof. Let m = [t/At\, then (4.1) provides the identity 



mAt (-L 



o J o 



(K u h)x -p(Qh)) (eh-^r) dxdt 



A '»V 



+1/2 



+ AtAx^^Up(/uV) i+1/2 



fc=0 i 



(t) 



(5.36) 



+ Ax^>A7 + 1 



Si "i t Si+i«i+i 



+1/2 



By virtue of the previous lemma, we know that the product terms involv- 
ing djAr^^-] converges to the product of the limits. Furthermore, from 
Lemmas 4.4 and 4.3, we know that the E™ and E™ terms converges to zero 
as h — > 0. Hence, the only term that might be problematic is the term 
involving Up(fwu). 

By adding and subtracting, we write 



AtAx£^U P (/S fc n fc ) i+1/2 



k=0 i 



(t) 



A* 



fc=0 i 
M 

AtAxYY,^^ u i+i/2 



k=0 i 
T r L 



J 



+ 



.M 



AtAx\YY.^)(ut 



+1/2 U i-1/2J 



fc=0 i 



where the last term is bounded as 

M 



/2 U i-l/2; 



k=0 i 



< ( A ^)IIS/i'"/i||L 2 (0,T;L2(o,L))ll(w/i)a;||L 2 (0,T;L 2 (0,L)) 

1+- — - 

< (Ax) 2 ^ll^h^H^^L^o.L))^ 

where the last inequality is a standard inverse estimate. The two previous 
equations tell us that 



lim AiAx Y Y Up ^ 



kZ7k„,k\ 
U U Ji+i/2 



fe=0 i 



JO 



gu dxdt. 



Hence, there is no problems with passing to the limit h — > in (5.36) and 
conclude the proof. □ 
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6. Strong convergence (proof of Theorem 2.4) 

In the previous section, we proved that the method almost converges to 
a weak solution of the compressible Navier-Stokes equations. To conclude 
convergence of the method, and thereby conclude the main theorem, it only 
remains to prove that 

p(g) = g a.e. 

This is the topic of this section. To make this identification, we shall prove 
that gh converges strongly a.e. Our strategy will be to adopt the continuous 
arguments of Lions [15] to the numerical setting. For this purpose, we shall 
need the following well-known lemma. The reader may consult [4] for a 
proof. 

Lemma 6.1. Let (g,u) satisfy the continuity equation (1.1) in the sense of 
distributions. If g G L 2 (0, T; L 2 (0, L)) and u G L 2 (0, T; Wq ,2 (0, L)), then 
for all B G C l {R+) such that B(g) G L 1 (0, T; L 1 (0, L)), 

B(g)t + (B(g)u) x + b(g)u x = 0, b(g) = gB'(g) - B(g), 

in the sense of distributions on [0,T) x [0,L]. 

From Lemma 5.3, we know that the limit (g, u) is a weak solution of 
the continuity equation and moreover that g G L 2 (0, T; L 2 (0, L)) and u G 
L 2 (0,T; Wq' 2 (0,L)). The previous lemma is then applicable. By setting 
B(z) = zlogz and integrating over the domain, we obtain 

/ glogg dx (t) = / g loggodx- / / gu x dx, (6.1) 
Jo Jo Jo Jo 

for any t G (0,T). 

Next, we set B(z) = zlogz in the renormalized continuity scheme (3.1), 
multiply with Ax, and sum over k = l,...,m, to obtain, for any m = 
1,...,M, 

rL pL pmAt pL 

/ Qh l °SQhdx< Q° h logQ° h dx- / g h (u h ) x dxdt, (6.2) 

Jo Jo Jo Jo 

where we have also used convexity of the map z i— >■ B(z) to conclude a sign 

on the error terms. 

By subtracting (6.1) from (6.2) and passing to the limit h — > 0, we deduce 

0<^y g log g - g log g dx^j (t) < J J gu x - gu^ dxdt, (6.3) 

for any t G (0,T). Here, the overbar denotes weak LMimit and the first 
inequality is a consequence of convexity. 

The remaining ingredient to prove strong convergence of the density is to 
prove that the right-hand side in (6.3) is negative: 

Lemma 6.2. Let (gh,Uh) be the numerical solution constructed using Defi- 
nition 2.1 and (2.3)-(2.4). Then, 

/ / gu x — gu x dxdt < 0. (6-4) 
Jo Jo 

As a consequence, Qh —> g a.e as h 0. 
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Proof. Let A n denote the Laplace operator with homogenous Neumann con- 
ditions on (0, L). From (5.1) and Lemma 3.5, we have that g S L°°(0, T; -L 7 (0, L)) 
and gt G L 2 (0,T;^-^(0,L)), respectively. As a consequence, 

M , 

g — ax, 



V = ^Ajy 1 



M~ 






-L 



L 



satisfies v x , vt £ L 2 (0, T; L 7 (0, L)). Hence, we can set v as test-function in 
(5.32) to obtain 

H r-L 



o J o 



fj,u x - p(g) 

t r-L 

o Jo 



g — — J dxat 



(gu)d t d x A N [g] + gu ( g 
L 



M 
L 



dxdt 



f L ( f L \ (6-5) 

- J dxA^loouoloo dx+ (J d x A^-[Qu]g dx\ (t) 

J J d t {d x A^ [gu]) g + gu 2 (p - ^-J dxdt 



where the last equality is integration by parts and the duality of the Neu- 
mann and Dirichlet Laplacians Ajy and Ad, respectively. 

To proceed, we note that cp = d x A D l [gu] is a valid test-function for the 
continuity equation (1.1). Since (g,u) is a weak solution of the continuity 
equation (Lemma 5.3), we deduce the identity 



/ / d t (d x A D l [gu\) g dxdt 
Jo Jo 



Jo 

t 



gud^Aplgu] dxdt 



L 

d x A D l [gouo]QQ dx 



t=o 



o Jo 



(gu) 2 dxdt 



d x A D [£> ^o]£o dx 



t=o 



By setting (6.6) in (6.5), we discover the identity 



o Jo 



\iu x - p(g) 



M 
L 

t 



o Jo 



) dxdt 

L 



gu dxdt + 



o 



(fj,u x -p(g)) g 
io Jo \ 

where the last equality is (5.35). 

Finally, we rewrite (6.7) in the form 



[iu x Q — fJ-UxQ dxdt 



M 
L 



dxAplQoUolQQ dx 

dxdt, 



t=o 



JO 



JO 



p{q)Q ~ p{q)q dxdt < 0, 



(6.6) 



(6.7) 



(6.8) 



where the last inequality follows from the convexity of z i— > p(z). This 
concludes the proof of (6.4). 
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By combining (6.8) and (6.3), we conclude that 

£>log g = glog q a.e on (0, T) x (0, L), 

and hence Qh —> Q a.e. □ 

Proof of Theorem 2.4: From Lemma 5.6, we have that (q,u) solves 

(gu) t + (gu 2 ) x = fj,u xx -p(q) x , 

in the sense of distributions. The previous lemma tell us that 

pjg) = p(g) a.e in (0, T) x (0, L). 

Hence, (g,u) is a weak solution of the compressible Navier-Stokes system 
(1.1)-(1.2). This concludes the proof of our main theorem (Theorem 2.4). 

□ 
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